library(tidyverse)
library(reshape2)
library(sf)
library(spdep)
library(tmap)
#library(OasisR) ## old package (OK for version R.2.1)
Import data for synthetic population located in ‘night’ cell
pop0, in ‘day’ cell pop1 and in ‘evening’ cell
pop2
pop_0 <- st_read( "data/population.gpkg",
layer="Population_0", quiet = TRUE)
pop_1 <- st_read("data/population.gpkg",
layer="Population_1", quiet = TRUE)
pop_2 <- st_read("data/population.gpkg",
layer="Population_2", quiet = TRUE)
Create new colums : gender (2 groups) ; Age (3 groups); Educ (3 Groups)
pop_0 <- pop_0 %>%
rowwise() %>%
mutate(
men = sum(across(starts_with('pop_1'))),
women = sum(across(starts_with('pop_2'))),
pop = men + women,
age1 = sum(across(starts_with(c('pop_1_1', 'pop_2_1')))),
age2 = sum(across(starts_with(c('pop_1_2', 'pop_2_2')))),
age3 = sum(across(starts_with(c('pop_1_3', 'pop_2_3')))),
educ1 = sum(across(ends_with('_1'))),
educ2 = sum(across(ends_with('_2'))),
educ3 = sum(across(ends_with('_3'))),
pct_educ3 = (educ3 / pop) * 100,
pct_age3 = (age3 / pop) * 100,
pct_women = (women / pop) * 100
)
pop_1 <- pop_1 %>%
rowwise() %>%
mutate(
men = sum(across(starts_with('pop_1'))),
women = sum(across(starts_with('pop_2'))),
age1 = sum(across(starts_with(c('pop_1_1', 'pop_2_1')))),
age2 = sum(across(starts_with(c('pop_1_2', 'pop_2_2')))),
age3 = sum(across(starts_with(c('pop_1_3', 'pop_2_3')))),
educ1 = sum(across(ends_with('_1'))),
educ2 = sum(across(ends_with('_2'))),
educ3 = sum(across(ends_with('_3')))
)
pop_2 <- pop_2 %>%
rowwise() %>%
mutate(
men = sum(across(starts_with('pop_1'))),
women = sum(across(starts_with('pop_2'))),
age1 = sum(across(starts_with(c('pop_1_1', 'pop_2_1')))),
age2 = sum(across(starts_with(c('pop_1_2', 'pop_2_2')))),
age3 = sum(across(starts_with(c('pop_1_3', 'pop_2_3')))),
educ1 = sum(across(ends_with('_1'))),
educ2 = sum(across(ends_with('_2'))),
educ3 = sum(across(ends_with('_3')))
)
Total Population density
bks_pop <- c(1, 10, 100, 1000, 10000, Inf)
tm_shape(pop_0) + tm_fill(col="pop", style="fixed",
breaks = bks_pop, palette="Greys") +
tm_legend(outside=TRUE)
Proportion of people with 2+ years of higher education
bks_edu <- c(0, 3.5, 7.7, 11.4, 17.2, 100)
tm_shape(pop_0) + tm_fill(col="pct_educ3", style="fixed",
breaks = bks_edu, palette="Reds") +
tm_legend(outside=TRUE)
Proportion of people aged 60+
bks_age <- c(0, 12.5, 17.1, 20.6, 25.1, 100)
tm_shape(pop_0) + tm_fill(col="pct_age3", style="fixed",
breaks = bks_age, palette="Blues") +
tm_legend(outside=TRUE)
Proportion of women
bks_sex <- c(0, 45.8, 50, 51.6, 54.5, 100)
tm_shape(pop_0) + tm_fill(col="pct_women", style="fixed",
breaks = bks_sex, palette="Greens") +
tm_legend(outside=TRUE)
Duncan’s segregation index is one-group form of dissimilarity index DI and measures the unevenness of a group distribution compared to the rest of the population. It can be interpreted as the share of the group that would have to move to achieve an even distribution compared to the rest of the population. You should not include a column with total population, because this will be interpreted as a group.
In night cells
dfpop_0 <- pop_0 %>% st_drop_geometry()
pop_0_Duncan_18categ<- as.data.frame(ISDuncan(dfpop_0[,2:19]))
colnames(pop_0_Duncan_18categ)[1] <- "Duncan_night"
rownames(pop_0_Duncan_18categ) <- paste (colnames(dfpop_0[ , 2:19]))
pop_0_Duncan_sex<- as.data.frame(ISDuncan(dfpop_0[,20:21]))
colnames(pop_0_Duncan_sex)[1] <- "Duncan_night"
rownames(pop_0_Duncan_sex) <- paste (colnames(dfpop_0[ , 20:21]))
pop_0_Duncan_age<- as.data.frame(ISDuncan(dfpop_0[,22:24]))
colnames(pop_0_Duncan_age)[1] <- "Duncan_night"
rownames(pop_0_Duncan_age) <- paste (colnames(dfpop_0[ , 22:24]))
pop_0_Duncan_educ<- as.data.frame(ISDuncan(dfpop_0[,25:27]))
colnames(pop_0_Duncan_educ)[1] <- "Duncan_night"
rownames(pop_0_Duncan_educ) <- paste (colnames(dfpop_0[ , 25:27]))
pop_0_all_Duncan <- bind_rows(pop_0_Duncan_18categ, pop_0_Duncan_sex, pop_0_Duncan_age, pop_0_Duncan_educ)
pop_0_all_Duncan_round <- pop_0_all_Duncan %>% mutate(across(starts_with("Duncan"), round, 3))
In day cells
dfpop_1 <- pop_1 %>% st_drop_geometry()
pop_1_Duncan_18categ<- as.data.frame(ISDuncan(dfpop_1[,2:19]))
colnames(pop_1_Duncan_18categ)[1] <- "Duncan_day"
rownames(pop_1_Duncan_18categ) <- paste (colnames(dfpop_1[ , 2:19]))
pop_1_Duncan_sex<- as.data.frame(ISDuncan(dfpop_1[,20:21]))
colnames(pop_1_Duncan_sex)[1] <- "Duncan_day"
rownames(pop_1_Duncan_sex) <- paste (colnames(dfpop_1[ , 20:21]))
pop_1_Duncan_age<- as.data.frame(ISDuncan(dfpop_1[,22:24]))
colnames(pop_1_Duncan_age)[1] <- "Duncan_day"
rownames(pop_1_Duncan_age) <- paste (colnames(dfpop_1[ , 22:24]))
pop_1_Duncan_educ<- as.data.frame(ISDuncan(dfpop_1[,25:27]))
colnames(pop_1_Duncan_educ)[1] <- "Duncan_day"
rownames(pop_1_Duncan_educ) <- paste (colnames(dfpop_1[ , 25:27]))
pop_1_all_Duncan <- bind_rows(pop_1_Duncan_18categ, pop_1_Duncan_sex, pop_1_Duncan_age, pop_1_Duncan_educ)
pop_1_all_Duncan_round <- pop_1_all_Duncan %>% mutate(across(starts_with("Duncan"), round, 3))
In evening cells
dfpop_2 <- pop_2 %>% st_drop_geometry()
pop_2_Duncan_18categ<- as.data.frame(ISDuncan(dfpop_2[,2:19]))
colnames(pop_2_Duncan_18categ)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_18categ) <- paste (colnames(dfpop_2[ , 2:19]))
pop_2_Duncan_sex<- as.data.frame(ISDuncan(dfpop_2[,20:21]))
colnames(pop_2_Duncan_sex)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_sex) <- paste (colnames(dfpop_2[ , 20:21]))
pop_2_Duncan_age<- as.data.frame(ISDuncan(dfpop_2[,22:24]))
colnames(pop_2_Duncan_age)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_age) <- paste (colnames(dfpop_2[ , 22:24]))
pop_2_Duncan_educ<- as.data.frame(ISDuncan(dfpop_2[,25:27]))
colnames(pop_2_Duncan_educ)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_educ) <- paste (colnames(dfpop_2[ , 25:27]))
pop_2_all_Duncan <- bind_rows(pop_2_Duncan_18categ, pop_2_Duncan_sex, pop_2_Duncan_age, pop_2_Duncan_educ)
pop_2_all_Duncan_round <- pop_2_all_Duncan %>% mutate(across(starts_with("Duncan"), round, 3))
Build table
pop_all_Duncan_round <- bind_cols(pop_0_all_Duncan_round, pop_1_all_Duncan_round, pop_2_all_Duncan_round)
pop_all_Duncan_round
In night cells
dfpop_0_prop <- dfpop_0/dfpop_0$pop
pop_0_geom <- pop_0[, -1:-28]
pop_0_prop <- bind_cols(pop_0_geom, dfpop_0_prop)
nb0 <- poly2nb(pop_0_prop, queen=TRUE)
lw0 <- nb2listw(nb0, style="W", zero.policy=TRUE)
list <- colnames(dfpop_0_prop[,2:27])
pop_0_moran <- data.frame(matrix(ncol = 2))
colnames(pop_0_moran) <- c("pop_night_IMoran" , "pop_night_Moranpvalue")
for (i in list)
{
pop_0_i <- pop_0_prop[,i]
colnames(pop_0_i)[1] ="ok"
pop_0_moran[i, 1] <- moran.test(pop_0_i$ok, lw0, zero.policy=TRUE, alternative="greater") [3]
pop_0_moran[i, 2] <- moran.test(pop_0_i$ok, lw0, zero.policy=TRUE, alternative="greater") [2]
}
pop_0_moran <- pop_0_moran[-1,] # Delete empty first line
pop_0_moran_round <- pop_0_moran %>% mutate(across(starts_with("pop"), round, 3))
pop_0_moran_round
In day cells
dfpop_1_prop <- dfpop_1/dfpop_1$pop
pop_1_geom <- pop_1[, -1:-28]
pop_1_prop <- bind_cols(pop_1_geom, dfpop_1_prop)
nb1 <- poly2nb(pop_1_prop, queen=TRUE)
lw1 <- nb2listw(nb1, style="W", zero.policy=TRUE)
list <- colnames(dfpop_1_prop[,2:27])
pop_1_moran <- data.frame(matrix(ncol = 2))
colnames(pop_1_moran) <- c( "pop_day_IMoran" , "pop_day_Moranpvalue")
for (i in list)
{
pop_1_i <- pop_1_prop[,i]
colnames(pop_1_i)[1] ="ok"
pop_1_moran[i, 1] <- moran.test(pop_1_i$ok, lw1, zero.policy=TRUE, alternative="greater") [3]
pop_1_moran[i, 2] <- moran.test(pop_1_i$ok, lw1, zero.policy=TRUE, alternative="greater") [2]
}
pop_1_moran <- pop_1_moran[-1,] # Delete empty first line
pop_1_moran_round <- pop_1_moran %>% mutate(across(starts_with("pop"), round, 3))
pop_1_moran_round
In evening cells
dfpop_2_prop <- dfpop_2/dfpop_2$pop
pop_2_geom <- pop_2[, -1:-28]
pop_2_prop <- bind_cols(pop_2_geom, dfpop_2_prop)
nb2 <- poly2nb(pop_2_prop, queen=TRUE)
lw2 <- nb2listw(nb2, style="W", zero.policy=TRUE)
list <- colnames(dfpop_2_prop[,2:27])
pop_2_moran <- data.frame(matrix(ncol = 2))
colnames(pop_2_moran) <- c( "pop_evening_IMoran" , "pop_evening_Moranpvalue")
for (i in list)
{
pop_2_i <- pop_2_prop[,i]
colnames(pop_2_i)[1] ="ok"
pop_2_moran[i, 1] <- moran.test(pop_2_i$ok, lw2, zero.policy=TRUE, alternative="greater") [3]
pop_2_moran[i, 2] <- moran.test(pop_2_i$ok, lw2, zero.policy=TRUE, alternative="greater") [2]
}
pop_2_moran <- pop_2_moran[-1,] # Delete empty first line
pop_2_moran_round <- pop_2_moran %>% mutate(across(starts_with("pop"), round, 3))
pop_2_moran_round
Combine Moran tables pop_0, pop_1, pop_2
pop_all_moran <- bind_cols(pop_0_moran_round, pop_1_moran_round, pop_2_moran_round)
data <- read.csv("data/sociodemo_distribution_healthy.csv", stringsAsFactors = F)
head(data)
## Sex Age Edu n_idf N_2002 N_2008 conso_5_2002 conso_5_2008
## 1 1 1 1 493871 52 148 0.01923077 0.02702703
## 2 1 1 2 496964 97 223 0.03092784 0.04932735
## 3 1 1 3 199476 84 137 0.03571429 0.07299270
## 4 1 2 1 1077670 241 355 0.04979253 0.08169014
## 5 1 2 2 662407 122 225 0.07377049 0.11111111
## 6 1 2 3 658465 118 180 0.05932203 0.16666667
n_idf = number of individuals of a particular group in
2012 in the Paris Region (Ile de France). N_2002 = number
of individuals of a particular group in the Health and Nutrition
Barometer survey in 2002. N_2008 = number of individuals of
a particular group in the Health and Nutrition Barometer in 2008.
conso_5_2002 = proportion of individuals eating 5+ portions
of fruit and vegetables a day in the Health and Nutrition Barometer in
2002. conso_5_2008 = proportion of individuals eating 5+
portions of fruit and vegetables a day in the Health and Nutrition
Barometer in 2008.
data$conso_5_2002 <- as.numeric(data$conso_5_2002)
data$conso_5_2008 <- as.numeric(data$conso_5_2008)
n_tot_idf <- sum(data$n_idf)
data$share_idf <- data$n_idf / n_tot_idf
data$h_idf_2002 <- round(data$n_idf * data$conso_5_2002,0)
data$h_idf_2008 <- round(data$n_idf * data$conso_5_2008,0)
data$h_baro_2002 <- round(data$N_2002 * data$conso_5_2002,0)
data$h_baro_2008 <- round(data$N_2008 * data$conso_5_2008,0)
total_prop <- function(data, year, sample){
if(sample == "idf"){
tp <-
sum(data[,paste("h", sample, year, sep="_")]) /
sum(data$n_idf)
} else {
tp <-
sum(data[,paste("h", sample, year, sep="_")]) /
sum(data[,paste("N", year, sep="_")])
}
return(tp)
}
Results with the Paris region as reference population
total_prop(data=data, year = 2002, sample = "idf")
## [1] 0.09568032
total_prop(data=data, year = 2008, sample = "idf")
## [1] 0.1205141
Results with the diet survey as reference population
total_prop(data=data, year = 2002, sample = "baro")
## [1] 0.1063931
total_prop(data=data, year = 2008, sample = "baro")
## [1] 0.1201574
n_tot_baro_2002 <- sum(data$N_2002)
n_tot_baro_2008 <- sum(data$N_2008)
EII <- function(data, year, sample){
if(sample == "idf"){
eii <- (
data[data$Sex == 1 & data$Age == 1 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 1 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 1, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 1 & data$Age == 2 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 2 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 2, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 1 & data$Age == 3 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 3 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 3, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 2 & data$Age == 1 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 1 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 1, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 2 & data$Age == 2 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 2 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 2, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 2 & data$Age == 3 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 3 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 3, "n_idf"] / n_tot_idf)
)
} else {
if(year == 2002){ n_tot_baro <- n_tot_baro_2002 }
if(year == 2008){ n_tot_baro <- n_tot_baro_2008 }
eii <- (
data[data$Sex == 1 & data$Age == 1 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 1 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 1, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 1 & data$Age == 2 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 2 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 2, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 1 & data$Age == 3 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 3 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 3, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 2 & data$Age == 1 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 1 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 1, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 2 & data$Age == 2 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 2 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 2, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 2 & data$Age == 3 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 3 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 3, paste0("N_", year)] / n_tot_baro)
)
}
return(eii)
}
Inequality results with the Paris region as reference population
EII(data=data, year = 2002, sample = "idf")
## [1] 1.407512
EII(data=data, year = 2008, sample = "idf")
## [1] 1.89476
Inequality results with the dietary survey as reference population
EII(data=data, year = 2002, sample = "baro")
## [1] 1.382566
EII(data=data, year = 2008, sample = "baro")
## [1] 1.905816
Importing calibration results
pareto <- read.csv2("data/nsga2_clean.csv")
head(pareto)
## evolution.generation evolution.evaluated evolution.samples maxProbaToSwitch
## 1 376800 376800 1 0.751
## 2 376800 376800 6 0.750
## 3 376800 376800 2 0.751
## 4 376800 376800 9 0.749
## 5 376800 376800 1 0.749
## 6 376800 376800 1 0.748
## constraintsStrength inertiaCoefficient healthyDietReward
## 1 0.100 1 0
## 2 0.100 1 0
## 3 0.098 1 0
## 4 0.097 1 0
## 5 0.097 1 0
## 6 0.097 1 0
## objective.deltaHealth objective.deltaSocialInequality X.Healthy_vs._2002
## 1 331,736,427 0.138 down
## 2 333,771,427 0.149 down
## 3 334,153,927 0.140 down
## 4 334,658,427 0.128 down
## 5 335,441,427 0.123 down
## 6 335,744,427 0.123 down
## X.Healthy_vs._2008 socialIneq_vs._2002 typo typo_name
## 1 down up down_down_up decrease%Healthy
## 2 down up down_down_up decrease%Healthy
## 3 down up down_down_up decrease%Healthy
## 4 down up down_down_up decrease%Healthy
## 5 down up down_down_up decrease%Healthy
## 6 down up down_down_up decrease%Healthy
## socialInequality
## 1 [1.7569191058549687]
## 2 [1.7525945924779105,1.7540192352458608,1.7427673698247963,1.749312240252607,1.7412857735082057,1.742133022717344]
## 3 [1.7625539632685843,1.745986394640198]
## 4 [1.7662218959212503,1.7696432627743401,1.7665457668386697,1.756434185222282,1.7667110959021646,1.7696451049204427,1.7694793140516023,1.763192441014371,1.7589183006338713]
## 5 [1.7713301959795948]
## 6 [1.7720880068092864]
## deltaSocialInequality
## 1 [0.13782823564503133]
## 2 [0.14215274902208952,0.14072810625413923,0.15197997167520372,0.14543510124739312,0.15346156799179433,0.15261431878265608]
## 3 [0.13219337823141575,0.14876094685980212]
## 4 [0.12852544557874968,0.12510407872565987,0.1282015746613303,0.138313156277718,0.12803624559783544,0.12510223657955732,0.1252680274483977,0.13155490048562912,0.13582904086612868]
## 5 [0.12341714552040517]
## 6 [0.12265933469071366]
## opinionMedian
## 1 [0.49988672137260437]
## 2 [0.4999333918094635,0.4999331533908844,0.49993330240249634,0.49993327260017395,0.4999333322048187,0.49993306398391724]
## 3 [0.4997502863407135,0.49975064396858215]
## 4 [0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5]
## 5 [0.5]
## 6 [0.499999463558197]
## numberOfHealthy
## 1 [834951]
## 2 [834948,832940,835222,834168,833731,834948]
## 3 [827118,826925]
## 4 [818522,817720,816927,815986,817833,817780,815267,817438,818191]
## 5 [817706]
## 6 [811529]
## deltaHealth
## 1 [331736.42711645283]
## 2 [332013.5408605285,334701.42711645283,332841.42711645283,332113.42711645283,335052.42711645283,334725.42711645283]
## 3 [333369.42711645283,334938.42711645283]
## 4 [333959.42711645283,333931.42711645283,335330.42711645283,334215.42711645283,336478.42711645283,336057.42711645283,334658.42711645283,335249.42711645283,334452.42711645283]
## 5 [335441.42711645283]
## 6 [335744.42711645283]
Drawing the Pareto front to select a parameter set (figure 4)
pareto$reliability.sample <- ifelse(pareto$evolution.samples>5, 2,1)
pareto$objective.deltaHealth <- as.numeric(gsub(",","",pareto$objective.deltaHealth))
front <- ggplot(data=pareto, aes(x = objective.deltaHealth, y=objective.deltaSocialInequality))
front +
geom_point(aes(group = typo_name, colour = typo_name,
size=reliability.sample)) +
scale_size_continuous(range = c(0.5,1.5)) +
scale_colour_manual(values = alpha(
c("#faa911", #orange
"#07b55e", #green
"#1088b0", #blue
"#aae860" #light green
), .5)) +
geom_vline(xintercept = 290000000) +
geom_hline(yintercept = 0.36) +
theme_bw()
Distribution of parameter values along the Pareto front (Appendix 6)
pareto$id <- rownames(pareto)
colnames(pareto)
## [1] "evolution.generation" "evolution.evaluated"
## [3] "evolution.samples" "maxProbaToSwitch"
## [5] "constraintsStrength" "inertiaCoefficient"
## [7] "healthyDietReward" "objective.deltaHealth"
## [9] "objective.deltaSocialInequality" "X.Healthy_vs._2002"
## [11] "X.Healthy_vs._2008" "socialIneq_vs._2002"
## [13] "typo" "typo_name"
## [15] "socialInequality" "deltaSocialInequality"
## [17] "opinionMedian" "numberOfHealthy"
## [19] "deltaHealth" "reliability.sample"
## [21] "id"
pareto_long <- melt(pareto[,c("maxProbaToSwitch","constraintsStrength" ,
"inertiaCoefficient","healthyDietReward",
"objective.deltaSocialInequality", "objective.deltaHealth",
"typo_name", "reliability.sample")],
id.vars=c("objective.deltaSocialInequality", "objective.deltaHealth",
"typo_name", "reliability.sample"))
ggplot(pareto_long, aes(objective.deltaHealth, value)) +
geom_point(aes(group = typo_name, colour = typo_name,
size=reliability.sample)) +
scale_size_continuous(range = c(0.5,1.5)) + scale_colour_manual(values = alpha(
c("#faa911", #orange
"#07b55e", #green
"#1088b0", #blue
"#aae860" #light green
), .5)) +
theme_bw() + theme(legend.position="none") +
facet_wrap(~variable, scales = "free_y", ncol = 1)
ggplot(pareto_long, aes(value, objective.deltaSocialInequality)) +
geom_point(aes(group = typo_name, colour = typo_name,
size=reliability.sample)) +
scale_size_continuous(range = c(0.5,1.5)) + scale_colour_manual(values = alpha(
c("#faa911", #orange
"#07b55e", #green
"#1088b0", #blue
"#aae860" #light green
), .5)) +
theme_bw() + theme(legend.position="none")+
facet_wrap(~variable, scales = "free_x", nrow = 1)
s0r <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg",
layer="Initial_0", quiet = TRUE)
s0o <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg",
layer="Initial_0", quiet = TRUE)
s1A <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg",
layer="Final", quiet = TRUE)
s1B <- st_read( "data/results_HigherProp_Scenario2_RandomPop_RandomMove_42.gpkg",
layer="Final", quiet = TRUE)
s2A <- st_read( "data/results_HigherProp_Scenario3_ObservedPop_NoMove_42.gpkg",
layer="Final", quiet = TRUE)
s2B <- st_read( "data/results_HigherProp_Scenario4_ObservedPop_RandomMove_42.gpkg",
layer="Final", quiet = TRUE)
s2C <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg",
layer="Final", quiet = TRUE)
bks <- c(0, 0.1, 0.125, 0.15, 0.175, Inf)
Initial Random
tm_shape(s0r) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Initial Observed
tm_shape(s0o) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 1A
tm_shape(s1A) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 1B
tm_shape(s1B) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 2A
tm_shape(s2A) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 2B
tm_shape(s2B) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 2C
tm_shape(s2C) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
duncan <- data.frame(matrix(ncol = 2, nrow=7))
colnames(duncan) <- c("type", "Duncan" )
dfs0r<- s0r %>% st_drop_geometry()
duncan[1,1] <- "s0r_night"
duncan[1,2]<- ISDuncan(dfs0r[,3:4]) [1]
dfs0o<- s0o %>% st_drop_geometry()
duncan[2,1] <- "s0o_night"
duncan[2,2]<- ISDuncan(dfs0o[,3:4]) [1]
dfs1A<- s1A %>% st_drop_geometry()
duncan[3,1] <- "s1A_final_night"
duncan[3,2]<- ISDuncan(dfs1A[,3:4]) [1]
dfs1B<- s1B %>% st_drop_geometry()
duncan[4,1] <- "s1B_final_night"
duncan[4,2]<- ISDuncan(dfs1B[,3:4]) [1]
dfs2A<- s2A %>% st_drop_geometry()
duncan[5,1] <- "s2A_final_night"
duncan[5,2]<- ISDuncan(dfs2A[,3:4]) [1]
dfs2B<- s2B %>% st_drop_geometry()
duncan[6,1] <- "s2B_final_night"
duncan[6,2]<- ISDuncan(dfs2B[,3:4]) [1]
dfs2C<- s2C %>% st_drop_geometry()
duncan[7,1] <- "s2C_final_night"
duncan[7,2]<- ISDuncan(dfs2C[,3:4]) [1]
duncan_round <- duncan %>% mutate(across(starts_with("Duncan"), round, 3))
duncan_round
moran <- data.frame(matrix(ncol = 3, nrow=7))
colnames(moran) <- c("type", "IMoran" , "pvalueMoran")
nbs0r <- poly2nb(s0r, queen=TRUE)
lws0r <- nb2listw(nbs0r, style="W", zero.policy=TRUE)
moran[1,1] <- "s0r_night"
moran[1,2] <- moran.test(s0r$X.propHealthy, lws0r, zero.policy=TRUE, alternative="greater") [3]
moran[1,3] <- moran.test(s0r$X.propHealthy, lws0r, zero.policy=TRUE, alternative="greater") [2]
nbs0o <- poly2nb(s0o, queen=TRUE)
lws0o <- nb2listw(nbs0o, style="W", zero.policy=TRUE)
moran[4,1] <- "s0o_night"
moran[4,2] <- moran.test(s0o$X.propHealthy, lws0o, zero.policy=TRUE, alternative="greater") [3]
moran[4,3] <- moran.test(s0o$X.propHealthy, lws0o, zero.policy=TRUE, alternative="greater") [2]
nbs1A <- poly2nb(s1A, queen=TRUE)
lws1A <- nb2listw(nbs1A, style="W", zero.policy=TRUE)
moran[2,1] <- "s1A_final_night"
moran[2,2] <- moran.test(s1A$X.propHealthy, lws1A, zero.policy=TRUE, alternative="greater") [3]
moran[2,3] <- moran.test(s1A$X.propHealthy, lws1A, zero.policy=TRUE, alternative="greater") [2]
nbs1B <- poly2nb(s1B, queen=TRUE)
lws1B <- nb2listw(nbs1B, style="W", zero.policy=TRUE)
moran[3,1] <- "s1B_final_night"
moran[3,2] <- moran.test(s1B$X.propHealthy, lws1B, zero.policy=TRUE, alternative="greater") [3]
moran[3,3] <- moran.test(s1B$X.propHealthy, lws1B, zero.policy=TRUE, alternative="greater") [2]
nbs2A <- poly2nb(s2A, queen=TRUE)
lws2A <- nb2listw(nbs2A, style="W", zero.policy=TRUE)
moran[5,1] <- "s2A_final_night"
moran[5,2] <- moran.test(s2A$X.propHealthy, lws2A, zero.policy=TRUE, alternative="greater") [3]
moran[5,3] <- moran.test(s2A$X.propHealthy, lws2A, zero.policy=TRUE, alternative="greater") [2]
nbs2B <- poly2nb(s2B, queen=TRUE)
lws2B <- nb2listw(nbs2B, style="W", zero.policy=TRUE)
moran[6,1] <- "s2B_final_night"
moran[6,2] <- moran.test(s2B$X.propHealthy, lws2B, zero.policy=TRUE, alternative="greater") [3]
moran[6,3] <- moran.test(s2B$X.propHealthy, lws2B, zero.policy=TRUE, alternative="greater") [2]
nbs2C <- poly2nb(s2C, queen=TRUE)
lws2C <- nb2listw(nbs2C, style="W", zero.policy=TRUE)
moran[7,1] <- "s2C_final_night"
moran[7,2] <- moran.test(s2C$X.propHealthy, lws2C, zero.policy=TRUE, alternative="greater") [3]
moran[7,3] <- moran.test(s2C$X.propHealthy, lws2C, zero.policy=TRUE, alternative="greater") [2]
moran_round <- moran %>% mutate(across(2:3, round, 3))
moran_round
## type IMoran pvalueMoran
## 1 s0r_night 0.010 0.061
## 2 s1A_final_night 0.010 0.057
## 3 s1B_final_night 0.016 0.008
## 4 s0o_night 0.006 0.175
## 5 s2A_final_night 0.033 0.000
## 6 s2B_final_night 0.024 0.000
## 7 s2C_final_night 0.009 0.091
s1A_day <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg", quiet = TRUE,
layer="Final_1")
s1B_day <- st_read( "data/results_HigherProp_Scenario2_RandomPop_RandomMove_42.gpkg", quiet = TRUE,
layer="Final_1")
s2A_day <- st_read( "data/results_HigherProp_Scenario3_ObservedPop_NoMove_42.gpkg", quiet = TRUE,
layer="Final_1")
s2B_day <- st_read( "data/results_HigherProp_Scenario4_ObservedPop_RandomMove_42.gpkg", quiet = TRUE,
layer="Final_1")
s2C_day <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg", quiet = TRUE,
layer="Final_1")
Duncan
duncan_day <- data.frame(matrix(ncol = 2, nrow=5))
colnames(duncan_day) <- c("type", "Duncan_dayCell" )
dfs1A_day<- s1A_day %>% st_drop_geometry()
duncan_day[1,1] <- "s1A_final_day"
duncan_day[1,2]<- ISDuncan(dfs1A_day[,3:4]) [1]
dfs1B_day<- s1B_day %>% st_drop_geometry()
duncan_day[2,1] <- "s1B_final_day"
duncan_day[2,2]<- ISDuncan(dfs1B_day[,3:4]) [1]
dfs2A_day<- s2A_day %>% st_drop_geometry()
duncan_day[3,1] <- "s2A_final_day"
duncan_day[3,2]<- ISDuncan(dfs2A_day[,3:4]) [1]
dfs2B_day<- s2B_day %>% st_drop_geometry()
duncan_day[4,1] <- "s2B_final_day"
duncan_day[4,2]<- ISDuncan(dfs2B_day[,3:4]) [1]
dfs2C_day<- s2C_day %>% st_drop_geometry()
duncan_day[5,1] <- "s2C_final_day"
duncan_day[5,2]<- ISDuncan(dfs2C_day[,3:4]) [1]
duncan_day <- duncan_day %>% mutate(across(starts_with("Duncan"), round, 3))
Moran
moran_day <- data.frame(matrix(ncol = 3, nrow=5))
colnames(moran_day) <- c("type", "IMoran_dayCell" , "pvalueMoran_dayCell")
nbs1A_day <- poly2nb(s1A_day, queen=TRUE)
lws1A_day <- nb2listw(nbs1A_day, style="W", zero.policy=TRUE)
moran_day[1,1] <- "s1A_final_day"
moran_day[1,2] <- moran.test(s1A_day$X.propHealthy, lws1A_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[1,3] <- moran.test(s1A_day$X.propHealthy, lws1A_day, zero.policy=TRUE, alternative="greater") [2]
nbs1B_day <- poly2nb(s1B_day, queen=TRUE)
lws1B_day <- nb2listw(nbs1B_day, style="W", zero.policy=TRUE)
moran_day[2,1] <- "s1B_final_day"
moran_day[2,2] <- moran.test(s1B_day$X.propHealthy, lws1B_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[2,3] <- moran.test(s1B_day$X.propHealthy, lws1B_day, zero.policy=TRUE, alternative="greater") [2]
nbs2A_day <- poly2nb(s2A_day, queen=TRUE)
lws2A_day <- nb2listw(nbs2A_day, style="W", zero.policy=TRUE)
moran_day[3,1] <- "s2A_final_day"
moran_day[3,2] <- moran.test(s2A_day$X.propHealthy, lws2A_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[3,3] <- moran.test(s2A_day$X.propHealthy, lws2A_day, zero.policy=TRUE, alternative="greater") [2]
nbs2B_day <- poly2nb(s2B_day, queen=TRUE)
lws2B_day <- nb2listw(nbs2B_day, style="W", zero.policy=TRUE)
moran_day[4,1] <- "s2B_final_day"
moran_day[4,2] <- moran.test(s2B_day$X.propHealthy, lws2B_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[4,3] <- moran.test(s2B_day$X.propHealthy, lws2B_day, zero.policy=TRUE, alternative="greater") [2]
nbs2C_day <- poly2nb(s2C_day, queen=TRUE)
lws2C_day <- nb2listw(nbs2C_day, style="W", zero.policy=TRUE)
moran_day[5,1] <- "s2C_final_day"
moran_day[5,2] <- moran.test(s2C_day$X.propHealthy, lws2C_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[5,3] <- moran.test(s2C_day$X.propHealthy, lws2C_day, zero.policy=TRUE, alternative="greater") [2]
moran_round_day <- moran_day %>% mutate(across(2:3, round, 3))
Combine Duncan and Moran files
Duncan_moran_dayCell <- left_join(duncan_day, moran_round_day, by = "type")
Duncan_moran_dayCell
s1A_evening <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg", quiet = TRUE,
layer="Final_2")
s1B_evening <- st_read( "data/results_HigherProp_Scenario2_RandomPop_RandomMove_42.gpkg", quiet = TRUE,
layer="Final_2")
s2A_evening <- st_read( "data/results_HigherProp_Scenario3_ObservedPop_NoMove_42.gpkg", quiet = TRUE,
layer="Final_2")
s2B_evening <- st_read( "data/results_HigherProp_Scenario4_ObservedPop_RandomMove_42.gpkg", quiet = TRUE,
layer="Final_2")
s2C_evening <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg", quiet = TRUE,
layer="Final_2")
Duncan
duncan_evening <- data.frame(matrix(ncol = 2, nrow=5))
colnames(duncan_evening) <- c("type", "Duncan_eveningCell" )
dfs1A_evening<- s1A_evening %>% st_drop_geometry()
duncan_evening[1,1] <- "s1A_final_evening"
duncan_evening[1,2]<- ISDuncan(dfs1A_evening[,3:4]) [1]
dfs1B_evening<- s1B_evening %>% st_drop_geometry()
duncan_evening[2,1] <- "s1B_final_evening"
duncan_evening[2,2]<- ISDuncan(dfs1B_evening[,3:4]) [1]
dfs2A_evening<- s2A_evening %>% st_drop_geometry()
duncan_evening[3,1] <- "s2A_final_evening"
duncan_evening[3,2]<- ISDuncan(dfs2A_evening[,3:4]) [1]
dfs2B_evening<- s2B_evening %>% st_drop_geometry()
duncan_evening[4,1] <- "s2B_final_evening"
duncan_evening[4,2]<- ISDuncan(dfs2B_evening[,3:4]) [1]
dfs2C_evening<- s2C_evening %>% st_drop_geometry()
duncan_evening[5,1] <- "s2C_final_evening"
duncan_evening[5,2]<- ISDuncan(dfs2C_evening[,3:4]) [1]
duncan_evening <- duncan_evening %>% mutate(across(starts_with("Duncan"), round, 3))
Moran
moran_evening <- data.frame(matrix(ncol = 3, nrow=5))
colnames(moran_evening) <- c("type", "IMoran_eveningCell" , "pvalueMoran_eveningCell")
nbs1A_evening <- poly2nb(s1A_evening, queen=TRUE)
lws1A_evening <- nb2listw(nbs1A_evening, style="W", zero.policy=TRUE)
moran_evening[1,1] <- "s1A_final_evening"
moran_evening[1,2] <- moran.test(s1A_evening$X.propHealthy, lws1A_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[1,3] <- moran.test(s1A_evening$X.propHealthy, lws1A_evening, zero.policy=TRUE, alternative="greater") [2]
nbs1B_evening <- poly2nb(s1B_evening, queen=TRUE)
lws1B_evening <- nb2listw(nbs1B_evening, style="W", zero.policy=TRUE)
moran_evening[2,1] <- "s1B_final_evening"
moran_evening[2,2] <- moran.test(s1B_evening$X.propHealthy, lws1B_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[2,3] <- moran.test(s1B_evening$X.propHealthy, lws1B_evening, zero.policy=TRUE, alternative="greater") [2]
nbs2A_evening <- poly2nb(s2A_evening, queen=TRUE)
lws2A_evening <- nb2listw(nbs2A_evening, style="W", zero.policy=TRUE)
moran_evening[3,1] <- "s2A_final_evening"
moran_evening[3,2] <- moran.test(s2A_evening$X.propHealthy, lws2A_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[3,3] <- moran.test(s2A_evening$X.propHealthy, lws2A_evening, zero.policy=TRUE, alternative="greater") [2]
nbs2B_evening <- poly2nb(s2B_evening, queen=TRUE)
lws2B_evening <- nb2listw(nbs2B_evening, style="W", zero.policy=TRUE)
moran_evening[4,1] <- "s2B_final_evening"
moran_evening[4,2] <- moran.test(s2B_evening$X.propHealthy, lws2B_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[4,3] <- moran.test(s2B_evening$X.propHealthy, lws2B_evening, zero.policy=TRUE, alternative="greater") [2]
nbs2C_evening <- poly2nb(s2C_evening, queen=TRUE)
lws2C_evening <- nb2listw(nbs2C_evening, style="W", zero.policy=TRUE)
moran_evening[5,1] <- "s2C_final_evening"
moran_evening[5,2] <- moran.test(s2C_evening$X.propHealthy, lws2C_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[5,3] <- moran.test(s2C_evening$X.propHealthy, lws2C_evening, zero.policy=TRUE, alternative="greater") [2]
moran_round_evening <- moran_evening %>% mutate(across(2:3, round, 3))
Combine Duncan and Moran files
Duncan_moran_eveningCell <- left_join(duncan_evening, moran_round_evening, by = "type")
Duncan_moran_eveningCell
library(cowplot)
library(tidyverse)
library(here)
phase <- read.csv("/Users/ccottineau/GitHub/Paper/Notebook_R_fig_tables/data/phase_48.csv")
long_table <- data.frame()
i <- 0
j <- 1
time <- data.frame(
day = rep("day",18),
nd = sort(rep(0:5,3)),
step = rep("step",18),
ns = rep(0:2,6)) |>
mutate(time = paste0(day, nd, "_", step, ns)) |>
select(time)
for(i in 1:nrow(phase)){
nsenario <- phase[i, "numOfScenario"]
EII <- t(
str_split(
str_replace(
str_replace( phase[i, "socialInequality"] , c("\\["), ""),
c("\\["), ""
),
",", simplify = T
)
)[1:18,]
healthy <- t(
str_split(
str_replace(
str_replace( phase[i, "numberOfHealthy"] , c("\\["), ""),
c("\\["), ""
),
",", simplify = T
)
)[1:18,]
long_table[j:(j+17),"time"] <- time
long_table[j:(j+17),"scenario"] <- rep(nsenario, 18)
long_table[j:(j+17),"socialInequality"] <- EII
long_table[j:(j+17),"numberOfHealthy"] <- healthy
j <- j+18
}
plot_table <- long_table |>
mutate(prop_healthy = as.numeric(numberOfHealthy) / 8768898,
EII = as.numeric(socialInequality),
scenario = as.factor(ifelse(scenario == 2, "1B",
ifelse(scenario == 5, "2C", NA))
)
)
ggplot(plot_table,
aes(x = time, y = EII, group = scenario,
colour = scenario)) +
geom_point()
specialRead <- function(x) {
read_csv(file=x, col_names = TRUE, col_types= "nnnnnnnnnn")
}
concatFiles <- function (base, folder){
fullPath = paste(base,folder,sep="")
list.files(path = fullPath, full.names = TRUE, pattern = "(\\d_\\d_\\d)") %>% lapply(specialRead) %>% bind_rows
}
recode <- function(x) {
step_to_day <- c (
'00' = 'day 0 step 1',
'01' = 'day 0 step 2',
'02' = 'day 0 step 3',
'10' = 'day 1 step 1',
'11' = 'day 1 step 2',
'12' = 'day 1 step 3',
'20' = 'day 2 step 1',
'21' = 'day 2 step 2',
'22' = 'day 2 step 3',
'30' = 'day 3 step 1',
'31' = 'day 3 step 2',
'32' = 'day 3 step 3',
'40' = 'day 4 step 1',
'41' = 'day 4 step 2',
'42' = 'day 4 step 3',
'50' = 'day 5 step 1',
'51' = 'day 5 step 2',
'52' = 'day 5 step 3'
)
return(step_to_day[x])
}
Simulated proportions of healthy agents, of social inequality (EII) and of rank-ordered inequalities (Et)
hp_sc5_health <- read_csv("data/HigherProp/health.csv")
hp_sc5_health$propHealthy <- hp_sc5_health$healthy / hp_sc5_health$effective
hpLong_sc5_health <- hp_sc5_health %>%
gather(key=facteur, value=valeur, propHealthy,avgOpinion,socialInequality,e ) %>%
mutate(ds= recode(paste0(day, slice)))
hpLong_sc5_health <- hpLong_sc5_health %>% mutate (type="hp")
both_sc5_health <- hpLong_sc5_health %>% mutate(type = as.factor(type))
p_hp_propHealthy <- both_sc5_health %>% filter(facteur=="propHealthy") %>% ggplot( aes(x=ds, y=valeur,colour=type) )+
theme_bw()+
theme(axis.text.x = element_text(angle = 45, vjust = 0.9, hjust=1.0)) +
theme(axis.text.x.bottom = element_text(margin=margin(t=0, r=0,b=10,l=0)) ) +
geom_point() +
theme(legend.position = "none") +
labs(y= "Value", x = "prop. of healthy")
p_hp_ineq <- both_sc5_health %>% filter(facteur=="socialInequality") %>% ggplot( aes(x=ds, y=valeur,colour=type)) +
theme_bw()+
theme(axis.text.x = element_text(angle = 45, vjust = 0.9, hjust=1.0)) +
theme(axis.text.x.bottom = element_text(margin=margin(t=0, r=0,b=10,l=0)) ) +
geom_point() +
theme(legend.position = "none") +
labs(y= "Value", x = "EducIneqIndex (EII)")
p_hp_e <- both_sc5_health %>% filter(facteur=="e") %>% ggplot( aes(x=ds, y=valeur,colour=type)) +
theme_bw()+
theme(axis.text.x = element_text(angle = 45, vjust = 0.9, hjust=1.0)) +
theme(axis.text.x.bottom = element_text(margin=margin(t=0, r=0,b=10,l=0)) ) +
geom_point() +
theme(legend.position = "none") +
labs(y= "Value", x = "Rank-ordered inequalities (Et)")
p_hp_all <- plot_grid(p_hp_propHealthy, p_hp_ineq, p_hp_e, labels = "AUTO",nrow=1, align="h")
p_hp_all
##ggsave("out/p_hp_all.pdf",p_hp_all, width=20)
Simulated proportion of healthy agents and simulated opinion across the 18 sociodemographic categories
sortie_sc5_higherprop <- concatFiles("data", "/HigherProp")
sortie_sc5_higherprop$propHealthy <- sortie_sc5_higherprop$healthy / sortie_sc5_higherprop$effective
sortieLong_sc5_hp <- gather(sortie_sc5_higherprop, facteur, valeur, c(propHealthy, avgOpinion) )
sortieLong_sc5_hp <- mutate(sortieLong_sc5_hp, "ds" = gsub(" ", "", paste(paste("[",paste(sortieLong_sc5_hp$day, sortieLong_sc5_hp$slice)),"]")))
sortieLong_sc5_hp <- mutate(sortieLong_sc5_hp, ds= recode(paste0(day, slice)))
age.l <- c("15-29 yrs", "30-59 yrs", "60-75 yrs")
names(age.l) <- c("1", "2", "3")
sex.l <- c("Male", "Female")
names(sex.l) <- c("1", "2")
p1 <- ggplot(data=sortieLong_sc5_hp %>% filter(facteur=="propHealthy"), aes(x=ds, y=valeur, group=educ, fill=educ, colour=as.character(educ)))
p1 <- p1 + scale_color_discrete(
name = "Education level",
labels = c("low", "middle", "high")
)
p1 <- p1 + theme_bw() + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
p1 <- p1 + labs(title="Proportion of healthy by gender, age and education level", y = "Value", x = "Time Slice")
p1 <- p1 + geom_line(size=0.5) + facet_grid( sex ~ age, labeller = labeller(sex = sex.l, age = age.l), scales = "fixed") #show.legend = FALSE
p1 <- p1 + coord_cartesian( ylim = c(0, 0.3))
p2 <- ggplot(data=sortieLong_sc5_hp %>% filter(facteur=="avgOpinion"), aes(x=ds, y=valeur, group=educ, fill=educ, colour=as.character(educ)))
p2 <- p2 + scale_color_discrete(
name = "Education level",
labels = c("low", "middle", "high")
)
p2 <- p2 + theme_bw() + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
p2 <- p2 + labs(title="Average Opinion by gender, age and education level", color = "Education level", y = "Value", x = "Time Slice")
p2 <- p2 + geom_line(size=0.5) + facet_grid(sex ~ age, labeller = labeller(sex = sex.l, age = age.l), scales = "fixed") #show.legend = FALSE
p2 <- p2 + coord_cartesian( ylim = c(0.3, 0.7))
propHealtyAndOpinion_higherprop <- plot_grid(
p1, p2,
labels = "AUTO", ncol = 1)
print(propHealtyAndOpinion_higherprop)
##ggsave("out/propHealtyAndOpinion_higherprop.pdf")
Social inequalities in healthy behaviour for the 5 scenarios (section 4.1)
Importing simulation results
Distribution of simulated values of EducIneqIndex (EII) for the five space-time scenarios (10 000 replications per scenario) (figure 5)
Overall proportions of agents with healthy behaviour for the 5 scenarios (section 4.2)
!! TO DO JULIE (Table 5)